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SUMMARY 

In this paper, we propose a numerical method for computing solutions to Biot’s fully dynamic model of 
incompressible saturated porous media fU. Our spatial discretization scheme is based on the three-field 
formulation (u-w-p) and the coupling of aTowest order Raviart-Thomas mixed element ^ for fluid variable 
fields {w, p ) and a nodal Galerkin finite element for skeleton variable field (u). These mixed spaces are 
constructed based on the natural topology of the variables; hence, are physically compatible and able to 
exactly model the kind of continuity which is expected. The method automatically satisfies the well known 
LBB (inf-sup) stability condition and avoids locking that usually occurs in the numerical computations in 
the incompressible limit and very low hydraulic conductivity. In contrast to the majority of approaches, 
our three-field formulation can fully capture dynamic behavior of porous media even in high frequency 
loading phenomena with considerable fluid acceleration such as liquefaction and biomechanics of porous 
tissues under rapid external loading. Moreover, we address the importance of consistent initial conditions 
for poroelasticity equations with the incompressibility constraint, which represent a system of differential 
algebraic equations. The energy balance equation is derived for the full porous medium and used to assess 
the stability and accuracy of our time integration. To highlight the capabilities of our method, a variety of 
numerical studies are provided including verification with analytical and boundary element solutions, wave 
propagation analyses, hydraulic conductivity effects on damping and frequency content, energy balance 
analyses, mass lumping considerations, effects of mesh pattern and size, and stability analyses. We also 
explain some discrepancies commonly found in dynamic poroelasticity results in the literature. Copyright 
© 2014 John Wiley & Sons, Ltd. 
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1. INTRODUCTION 


Dynamics of porous media are of interest in a variety of fields including mechanics of saturated 
and partially saturated soil, porous biological materials and petroleum reservoirs. The dynamic 
response of saturated soil under earthquake and the subsequent soil liquefaction process is subject 
of many studies m- Additionally, the ability to accurately model the behavior of rock/soil 
and fluids in petroleum reservoirs is essential for maximizing oil recovery and assessing the 
performance of processes such as wellbore stability and subsidence. For this also, coupled modeling 
of geomechanics and fluid flow is required @0- Similarly in biomechanics, mechanical response 
of hydrated soft tissues such as brain, skin, and articular cartilage |[8f[T^, bones p3| , and even 
individual cells 1 1^ 1^ involves coupled response of a skeleton and pore fluid. Furthermore, all 
the results of dynamic poromechanics can be applied to dynamic thermomechanics based on the 
mathematical analogy between the formulations 
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In this paper, we propose a numerical method for dynamic poroelasticity, starting from Biot’s 
formulation Q. The porous medium is considered as a solid skeleton saturated with fluid. While 
the skeleton itself is deformable and elastic, the material that makes up the skeleton and the fluid 
are modeled as incompressible. This is reasonable for many of the physical problems described 
above. We also restrict the formulation to linearized kinematics. Alternate formulations based on the 
Theory of Porous Media p8| have been shown to be equivalent in this realm of incompressible solid 
and fluid materials and linear kinematics 1191. Our focus is on fully dynamic problems, i.e., inertia 
effects in both the fluid and the skeleton have been taken into account. We note that this is in contrast 
to most previous work, where fluid acceleration terms are fully or partially neglected 


11 201 . 


Some of the specific features of our approach are as follows: (a) In contrast to much of the 
literature, we use a full three-field formulation in which the skeleton displacement (rt), pore pressure 
(p) and pore-fluid relative velocity (w) fields are all approximated independently and calculated 
directly. As a result (see section ED, the formulation is applicable to problems with considerable 
fluid acceleration, for instance in liquefaction analysis 0 and biomechanics of porous tissues under 
rapid external loading p0l[T^ . Moreover, greater accuracy is attained compared to formulations 
where some fields are computed by post processing techniques as reported in | |^[^ . (b) We 
employ a topology-motivated spatial discretization based on the lowest order Raviart-Thomas (RTq) 
mixed element Q to approximate the flow variables and a linear/quadratic standard Galerkin finite 
element for the skeleton displacement. This mixed approach automatically enforces the precise 
kind of continuity (local mass conservation) that is required. Additionally, in the incompressible 
limit (when material is modeled as incompressible and under low hydraulic conductivity) and the 
rigid skeleton limit, the proposed method satisfies the well known LBB (inf-sup) condition | |2^[^ , 
avoiding associated locking effects and the need for any additional stabilization techniques, (c) 
Taking a differential algebraic equations (DAE) point of view, we demonstrate the importance of 
consistent initial conditions for poroelasticity equations with the incompressibility constraint, (d) 
The energy balance equation for the whole porous medium is also derived and used to confirm 
the stability and accuracy of the Newmark time integration method with particular emphasis on 
consistent initial condition. 

In many problems, such as soil mechanics and biomechanics of soft tissues, the fluid and skeleton 
material could be reasonably modeled as being incompressible. However when developing a finite 
element formulation, the incompressibility condition imposes a certain constraint (see section 3.2 1 , 
which in turn gives rise to stability concerns. The discrete form of the Ladyzenskaja-Babuska- 
Brezzi (LBB) condition provides the necessary and sufficient criterion for a stable finite element 
formulation. To satisfy this condition, several approaches including stabilized methods l|6 25 -301, 
discontinuous Galerkin 1M}|37), special elements such as Taylor-Hood elements p8) or MINI 
elements p^ , and mixed finite elements @E3ED have been proposed in the literature. As we 
discuss in section [T2] in a finite element formulation for poroelasticity, the LBB condition must be 

1. In this paper. 


satisfied in two limiting cases — low hydraulic conductivity and rigid skeleton | 
we develop a finite element formulation using a two fold approach to satisfy the LBB condition — 
(i) We discretize the pressure (p) and velocity (w) fields of the pore-fluid flow using Raviart-Thomas 
(RT) elements Q. These elements are topology-inspired, associating velocity degrees of freedom 
(DOL) with element edges and pressure DOL with the element interior, and result in automatically 
satisfying the LBB condition in the rigid skeleton limit, (ii) We discretize the skeleton displacement 
field (u) using standard Galerkin finite elements in such a way that together with the fluid pressure 
field, the LBB condition in the low hydraulic conductivity limit is satisfied in the same manner as in 
incompressible elasticity problems p^[4^ . In the context of RT elements, we note other topology- 
inspired finite elements p9] 44 1E3- 

Most finite elements of the three-field category for porous media are obtained using standard 
nodal Galerkin discretization of all three fields pl |T0l[^[46H5^ . However, these finite elements 
usually fail to fulfill the LBB condition. As described above, in this paper, we exploit a mixed 
finite element approach to deal with the challenges that conventional node-based Galerkin methods 
are not able to handle easily. The main advantage of mixed elements is their ability to exactly 
model the kind of continuity which is expected i.e., continuity of the normal velocity across element 
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edges which results in the desirable element-wise (local) mass conservation. Mixed elements form a 
structure that imitates, at the discrete level, all the operators that exist at the continuous level PDEs 
(such as divergence and gradient). In addition, boundary conditions can be applied more naturally. 
For these reasons, mixed elements such as RT elements Q, Nedelec elements m or a variety of 
others summarized in 1411 are widely used in solid and fluid mechanics problems of incompressible 
elasticity and fluid flow. However, their application in coupled problems (simultaneous solution) has 
been limited. Phillips and Wheeler | ^ 541 and Lipnikov | |55| developed a formulation based on a 
combination of the lowest order Raviart-Thomas (RTg) mixed finite element for the w and p fields 
and the standard nodal linear Galerkin finite element for the u field for quasistatic poroelasticity 
problems. We call this element Pi — RTg. They elaborated on the solvability and error estimates of 
this finite element for quasistatic consolidation problems. In this paper, we extend the application 
of this Pi — RTg element to dynamic problems. In doing so, we are able to effectively capture wave 
propagation effects (section [6T] l. We also observe that careful consideration of the finite element 
mesh is required to obtain the proper relationship between hydraulic conductivity and the damping 
of the mechanical response (section [6^ . This relationship has not been captured correctly by many 
past results in the literature. Furthermore, we find that while stability in terms of the FBB condition 
is achieved for typical ranges of hydraulic conductivity, in the limiting case of very small hydraulic 
conductivity, the Pi — RTg element exhibits instability in the form of checkerboard pattern in the p 
field. For reasons described in section [6!3.2[ this effect (also referred to as locking in this context), 
is more pronounced in dynamic problems than previously observed in quasistatic problems | |5^ . 
To resolve this issue, we propose a P 2 — RTg element which employs a quadratic Galerkin finite 
element to estimate the space of u. The element is shown to pass the numerical inf-sup test and 
therefore the stability criterion (see section 6.3.2 1 . 

This paper is organized as follows. First we summarize Biot’s mathematical model for porous 
media in section In section we motivate and present our mixed finite element formulation; 
show how the poroelasticity equations reduce to two classical models in the limiting case, and 
discuss the associated stability considerations; we develop our spatial discretization scheme based 
upon the three-field formulation (u-w-p) and describe the proposed finite element formulation using 
RT elements for the w-p fields and nodal Galerikin element for the u field. In section|^ we discuss 
time integration, in particular addressing the need for consistent initial conditions; we also present 
the energy balance equation for the full porous medium, which we later use to assess the stability 
and accuracy of time integration. We present the details of implementation in section including 
aspects involving unconventional edge degrees of freedom. To demonstrate the performance of this 
approach, we consider detailed numerical examples in section These include verification with 
analytical and boundary element solutions, wave propagation analysis, energy balance analyses, 
mass lumping considerations, effects of mesh pattern and size, hydraulic conductivity effects on 
damping and frequency, spurious pressure modes and locking, and stability analysis by means of 
inf-sup tests. We also comment on some discrepancies found in the results in the literature studies. 


2. GOVERNING EQUATIONS 

We describe Biot’s equations Q under the following conditions: (i) fluid and skeleton materials are 
incompressible; (ii) there is a single fluid phase; (iii) skeleton has a linear elastic behavior; (iv) a 
linear form of Darcy’s law (i.e,. constant hydraulic conductivity) is used; (v) a linearized kinematics 
formulation is applied which also results in constant porosity; (vi) the additional apparent mas^ 
used in Biot’s equations is neglected; (vii) to focus on the essential features of our formulation, 
we do not consider fluid sources or sinks, body forces or gravity effects, inclusion of which is 


talso called the added mass; in porous medium, when the fluid (skeleton) is moving, a force must be exerted on the 
skeleton (fluid) to prevent an average displacement of the latter which is defined through a mass coupling coefficient. 
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(b) 


Figure 1. (a) Saturated porous medium with primary variables, (b) Domain D corresponds to the skeleton 
boundaries is shown on the left while at the same time, D corresponds to the fluid boundaries is shown on 
the right. In each figure boundary F is divided into two disjoint subsets which cover the entire boundary of 

the domain; i.e., F^i U Fjv = F and Fp U Fq = F. 


straightforward. Accordingly, the linearized governing equations read [ST) 

pil + PfW — V ■ {a — pI) =0 

Pfu-\ — w + -—w + y[p) = u ( 1 ) 

nf iFh 

V • (m) + V • {w) =0 

These represent the momentum balance of the porous medium, the dynamic extension of Darcy’s 
law, and the mass balance equation (continuity equation) respectively. The three basic fields are 
skeleton displacement, u, pore fluid pressure, p, and fluid Darcy velocity, w, as shown in Figure 
[Ta| Darcy velocity is defined as the relative rate of discharge per unit area of the medium, 
w = rifjwf — ii), where uf refers to the fluid absolute velocity. V-, V, and dot denote the gradient 
operator, the divergence operator, and the derivative with respect to time, t, respectively, pf is the 
fluid density, while p = rif pf + (1 — nf)ps is the average density of the medium, with ps and rtf being 
solid phase density and porosity respectively. I, g, and ATh represent the identity matrix, acceleration 
due to gravity, and the hydraulic conductivity (its associated matrix is assumed to be isotropic). 

The concept of effective stress is used as a principle that controls the constitutive behavior of a 
porous medium: cr‘ = cr — ap I where the total stress is composed of effective stress of skeleton 
cr and the pore pressure p. Here, a denotes the Biot-Willis coefficient which is commonly taken 
as 1 for soil material. For an isotropic elastic skeleton, the effective stress is given by the usual 
constitutive equation 

a = 2Ge + XI trace(e) (2) 

where A and G are the Lame’s first and second constant (shear modulus) of the skeleton under 
drained condition, respectively, e = |(Vm + Vrt^) is the strain. 

Equation 0 is a system of coupled partial differential equations with u, w, and p as the three 
unknown fields on the domain D. The following are appropriate boundary conditions 

u = Mprsc over Fd and a'^ ■ fi = Tprsc over F m (skeleton boundary condition) 
w ■ h = gprsc over F^ and p = pprsc over Tp (fluid boundary condition) 

where Uprsc, Tprsc,Pprsc ^ud gprsc are the prescribed skeleton displacement, traction, pressure and 
volumetric flux, and n is the outward unit normal to the boundary F. We note that for each of 
the flow and skeleton, the boundary is decomposed into two disjoint subsets which cover the entire 
boundary F as shown in Figure [Tbl 
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3. FINITE ELEMENT EORMULATION 

The coupled nature of the poroelastic governing equations 0 precludes analytical solutions, except 
in some special cases like those presented in the seminal paper of Rice and Cleary | |58) . In 
most cases, a numerical approximation must be computed. Eor this purpose, typically a finite 
element approach is used for spatial discretization, resulting in a system of differential algebraic 
equations (DAE) that can be solved by a time integration scheme. In this section, we discuss spatial 
discretization, specifically our approach of combining RT elements for the w-p fields with standard 
nodal Galerkin elements for the u field. Eirst we motivate our choice of a three-field formulation. 
Then we discuss the stability concerns that arise. Einally, we present the finite element formulation. 


3.1. Three-field formulation 

Zienkiewicz and Shiomi 0 first discussed various possible finite element formulations for Biot’s 
dynamic poroelastcity model with the notion of incompressibility. These can be classified into three- 
field (u-w-p) formulations and two-field {u-p or u-w) formulations. Zienkiewicz and Shiomi 0 did 
not present an implementation of a three-field formulation. A majority of current implementations 
are based on two-field formulations obtained by eliminating either the fluid Darcy velocity w, 
resulting in a u-p formulation ||5||II||20||37||59f[62| or the pressure, leading to a u-w formulation 
0|9 63-661. These methods are efficient in terms of the computational cost as they reduce the 
number of DDEs. However, we pursue here a three-field formulation motivated by the following 
considerations. 

Generality - When the fluid and the solid material making up the skeleton are modeled as 
incompressible, the pressure term does not appear in the mass conservation equation 03. 
Consequently, the fluid Darcy velocity, w, cannot be directly eliminated from the Darcy’s law 02 
02g. Therefore, to obtain a u-p formulation, either (a) penalty procedures must be used in the mass 
conservation equation |3p3-65), or (b) the fluid relative acceleration, w, (and hence inertia) must be 
neglected. The latter is a simplification, first introduced by Zienkiewicz and coworkers 0|59), has 
been the most common means of obtaining a u-p formulation |[5] [TT][^[T7][59H6^ . This approach 
has been found to be typically valid for the slow and medium frequency level phenomena |59 [60) . 
However, in a variety of problems, it is desirable to take into account all acceleration terms; which 
is only possible by exploiting the three-field formulation, without using a penalty formulation. Eor 
instance using verification and validation procedures, Jeremic and coworkers 0 ^ 671 show that 
a three-field formulation is required to realistically model the saturated soil-foundation-structure 
interaction in earthquake as well as the liquefaction and cyclic mobility phenomena. They conclude 
that a full-field method can appropriately model the physical damping in the dynamic poroelasticity 
problem during seismic events by directly taking into account the fluid and skeleton interaction 
through Darcy’s equation. Accordingly, we could model and compute the velocity proportional 
damping energy as shown in equation ( |2^ whereas in other studies an artificial damping has to be 
used 0. Another area of interest is soft biological tissues that are often modeled mechanically as 
poroelastic media, and experience very rapid external forces m This situation creates significant 
inertia forces in the fluid phase that is not negligible and requires use of a full formulation |[^. 
An improved performance of the full-field mixed finite element formulation over the two-field 
formulation is also reported in pO|[TT| for large deformation simulations and in | |68) for modeling 
of bones subjected to Osteoporosis disease. Presence of the fluid velocity field in an explicit fashion 
will also allow the formulation to be more readily extended for coupled transport phenomena. 
Accuracy - The u-p formulation amounts to solving for pressure and skeleton displacement as 
primary variable fields. Typically post-processing techniques are required to compute the velocity 
field from the gradient of the pressure field, and is therefore of lower accuracy | 2Tp2 |. Eurthermore, 
both pressure and flux boundary conditions can be applied directly with a three-field formulation, 
but require special treatment if a two-field formulation is used 
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3.2. Nature of stability considerations 

To understand the nature of stability issues that are relevant in a three-field finite element 
formulation, we describe two limiting conditions in terms of hydraulic conductivity and rigidity of 
skeleton similar to ||6 421. Both of these limiting cases give rise to saddle point problems. Without 
loss of generality, for these limiting cases we only consider the static version of the poroelasticity 
problem where all the inertia terms are zercj^ 

Case (i) - Low hydraulic conductivity: In the impermeable limit (iTh —> 0), the static version of 
Q reduces to 


2GV^u-Vp = 0 
V • u =0 


(4) 


where is the Laplacian operator. Note V • ft = 0 is replaced by V • m = 0 as they are equivalent 
if the initial conditions satisfy. We recognize that these are the equations of the classical 
incompressible elasticity problem. Thus in the impermeable limit, there is no flow and the skeleton 
becomes incompressible but deformable. 

Case (ii) - Rigid skeleton: The other limiting case is when the skeleton matrix become rigid 
(which implies rt —> 0) and results in 


w + — Vp = 0 
Pf9 


(5) 


V • w 


= 0 


These are the classical Darcy’s flow equations. We note that this second limiting case applies only 
to a three-field formulation. 

Both these limit cases are saddle-point problems fTO) . The coefficient matrix associated with 
the discrete form of the saddle point problem has a zero-block structure on lower diagonal (see 
equation ( |^ for example) that creates stability challenges. There are two stability criteria that 
provide the necessary and sufficient condition for well-posedness of such problems, namely the 
ellipticity condition and the LBB condition (also called inf-sup condition). We refer to |70 IIZD for 
a thorough description on this subject. The discrete version of the ellipticity condition is satisfied 
for any choice of finite element with positive definite mass and stiffness matrices upon application 
of boundary conditions. The LBB stability condition is more important and challenging. Often, 
standard finite element approaches presented in the literature fail to satisfy this condition, leading 
to instability manifested as checkerboard patterns in the pressure field poor convergence behavior. 


We provide a comparison of stable and unstable finite element results in section 6.3 

Since the poroelasticity model degenerates to the saddle point problems Q and (0, a finite 
element formulation must satisfy the stability conditions associated with both these limiting cases. 
Thus the poroelasticity problem has the additional requirement that both the u-p pair and the w-p 


pair must be able to satisfy the stability conditions simultaneously. We show in section 6.3.1 that 
the proposed P 2 — RTq element qualifies. 

3.3. Finite element discretization 

We start from the weak-form representation of 0, 

p / 5u^ildVL + pf Su^wdil + / {N6u)^{< 7 — pl)d^l — / Su^ TprscdT + 

J Vt J Q, J Q, 

iw^ildQ + — [ Siv^w dfl + / Sw^wdil — [ ■ n)p„rsc dT+ 

Jn dih Jn Jn Jr^ 

/ V • (u) dQ + / Sp^ ■ {w) dn = 0 
Jfi Jq. 


Pf 


( 6 ) 


^The dynamic version affects only on ellipticity condition. 
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where 5u, 5w, and 6p are the skeleton displacement, fluid relative velocity, and pressure test 
functions respectively. We recognize now that the functions u, w, and p must belong to the spaces 

U = {u & : u = Upi-sc over Tu} 

W = {w G i7(div; Q) : w ■ fi = gprsc over r^} (7) 

v = {pGL^{n)} 


respectively. Here is the space of square integrable functions on Q, is the usual 

Sobolev space in which each function and its first order derivative belong to and i7(div; H) 

represents the space of vector-valued functions that belong to L^(H) and whose divergences belong 
to L^(H) as well. Further elaboration on these spaces can be found in Jti) and also in ||2 53 
Similarly, the test functions of Su, 6w, and Sp belong to the same spaces as their corresponding 
variable fields, except that instead of Uprsc and gp^c a zero is placed. We call these test spaces Uq, 
Wo, and Vq. 


As discussed above, with the goal of satisfying the appropriate stability conditions, we 
simultaneously apply the lowest order Raviart-Thomas mixed finite element (RTq) Q for flow 
variable fields and the standard linear/quadratic Galerkin finite elements for skeleton variable fields. 
Let the domain H be partitioned with T = {71, • • • ,Tnci} into N^i non-overlapping triangular 
elements. We denote by and N^g number of nodes and edges. Also position vector is a; = 

[a; y]^ ■ 

Skeleton displacement discretization: We discretize the skeleton displacement field using 
standard linear/quadratic shape functions over triangles, Ni{x), 


Nn 

u{x,t) = '^^Ni{x)ui{t) = N{x)u{t) (8) 

i=l 


We note that the components of vector u{t) are skeleton displacement DOFs, Ui, that are 
displacements at the element nodes. The skeleton and fluid DOFs are illustrated in Figure We 
call the combination of linear Galerkin finite element with RTq finite element. Pi — RTq, and the 
combination of quadratic Galerkin finite element with RTq finite element, P 2 — RTq. We propose 
the latter finite element to resolve the stability issues that the linear element creates in the limiting 


case of very small hydraulic conductivity as elaborated in section 6.3 

Fluid velocity discretization: 


w(x,t) ^'^Wj{x)qj{t) = W{x)q{t), Wj{x) 
i=i 


^ if X G 'Trn 
0 otherwise 


(9) 


Based on RTq element, Darcy velocity DOFs qj (Figure]^ are values of the normal velocity at the 
midpoints of edges of element m (in fact, normal velocity is constant along each edge), defined as 
ijj = Wj ■ fimj where Wj is fluid velocity vector defined on edge j. Also fimj denotes the outward 
unit normal vector to edge j of element m, while is the unit normal vector to edge j chosen with a 
global fixed orientation. This global convention is such that Cj points from element Tm with smaller 
index into its adjacent element with larger index (see Figurej^. Therefor, based on this convention, 
edges on the boundaries always have an outward global unit normal. In equation ( fT0| i Xj is the 
position vector of the vertex node opposite the j* edge of element m and Smj = nmj ■ Sj is a sign 
convention equal to -fl if ij points outward and otherwise is —1. The effect of this sign convention 
can be seen in Figur^^ where the velocity interpolation function, Wj{x), for the common edge (j) 
of the two adjacent elements is illustrated. Due to the convention, the element with smaller index Tm 
has Smj = 1 and the one with larger index Tm+i has Sm+ij = —1- Moreover, shown in Figure]^ is 
the normal component of velocity along the common edge of two adjacent elements illustrating the 
continuity of such normal component. Also Am and Ij refer to the area of element m and the length 
of edge j, respectively. 
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Figure 2. Fluid and skeleton DOFs in the coupled Galerkin/RTg elements - (a) linear Galerkin/RTo 
element, referred as Pi — RTq, (b) quadratic Galerkin/RTo element referred as P 2 — RTg. Both elements 
are constructed on the natural topology of the variable helds, i.e., skeleton displacement fields are defined 
on the nodes using a standard Galerkin approximation, while fluid variable fields are defined on the element 
edges (velocity) and volumes (pressure) using RTq element. We note that Darcys velocity DOFs are the 
values of the normal velocity at the midpoints of edges of elements and the pressure DOFs are set to be the 

constant elemental pressure. 



Figure 3. (a) Two enumerated adjacent elements Tm and T^+i that share the same edge j. In this figure, 
hmj and fim+ij are the outward unit normal vector to edge j of elements m and m + 1, respectively, while 
ij is the unit normal vector to edge j chosen with a global fixed orientation. This is such that Cj is equal to 
the outward unit normal vector of element with smaller index, i.e., Tm (and hence points into Tm+i). Also 
based on this convention, the element with smaller index Tm has Smj ~ 1 whereas the one with larger index 
Tm+i has Sm+ij = —1- (b) velocity interpolation function for the common edge of two adjacent elements, 
Wj{x). (c) normal velocity along the common edge of two adjacent elements. Note the continuity of this 
component of velocity across element edge which results in element-wise mass conservation. 


Pressure discretization: Pressure DOF pm is set to be constant elemental pressure, thus its 
interpolation function Hm (a:) is a polynomial of degree 0 that takes a constant value of 1 at element 
m and 0 at all other elements. 


Ael 

p{x,t) = ^ Hm{x)pm{t) = H{x)p{t), 

m—l 




if X GTm 
otherwise 


( 10 ) 


Fluid Darcy velocity and pressure DOFs are collected in vector q{t) and p{t), respectively and are 
shown in Figure]^ 
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Given the described finite element approximation, we formulate a discrete version of the weak 
problem (|^ by substituting approximation functions of u, in it. 


/ 

Jn 


,/ 

Jq 


/ 

Ja 


Su'N'Nudn-\-pf Su' Wqdn + / Su'B'CBudQ- / Su' B' 6HpdQ = 


j 

Ja 


L 


5 u' N' 


Pf 


/ dq^W^Nudn 

Ja 


j 

Ja 


rif 

L 


[ Sq^W^Wqdn 
Jn 


Pi 9 
A'h 


/ 

Ja 


L 


6 p^Budn + / 5 p^{W ■W)qdVi = 0 


Sq^W^Wqdn- / 6q^{SJ-W)^Hpdn = - 6q^{W^ ■h)p^„^dr 


( 11 ) 


The solution of the discretized problem must satisfies the integral equations for all test functions 
((Jn, 6w, 6p) gUq X VVq x Vq. In equation (ED 6 denotes the Kronecker delta in the vector format, 
B{x) is the linearized deformation-displacement matrix, and C is the elastic constitutive matrix under 
the drained condition. Note that we dropped the dependence of functions to x and t due to clarity. 
The test functions and the shape functions coincide in this formulations. Moreover, equation ( [TT] ) 
represent a consistent variational formulation in which the same shape functions are applied for 
higher time derivatives. An alternative formulation is discussed in section 6.1.5 The prescribed 
skeleton displacement and fluid flux are imposed in a strong way, i.e., as essential boundary 
conditions, while the prescribed traction (Tprsc) and pressure (pprsc) are entered the variational 
form as natural boundary conditions. Also, the constitutive equation (|^ is used in strong form 
as substitute in ([T]li. The compact form of ED is the following semi-discrete expression 


Q^ii +B^q 

where the following terms are defined 


Mil +Mfq+Ku —Qp = F 

Mf^il + Aq +^^Aq — Bp = ¥ 
Kh 

= 0 


M = p f N' Ndn 
Jn 


A=^ W^Wdn 


M{ = pf / N^WdFl K= / B^CBdn 


I 

Jn 


Q= / B'SHdn 
Jn 


B= / {\/ -W)' HdFl 
Jn 


( 12 ) 


(13) 


rp,se dr ¥ = - ■ n)pprsc dr 


We reiterate that such topology-based finite element discretization is physically compatible, and as 
a result all the above constant matrices (at the discrete level) imitates their corresponding operators 
at the continuous level in ([TJ; for instance B is the discrete version of divergence operator whose 
m}'^ column represents the m* element where all entries are zeros except the rows representing the 
free edges belongs to that element. In those entries depending on the flux convention, ±1 is placed. 
Further elaboration on building these matrices is provided in Section Also using the described 
finite element, specification of pressure and flux boundary conditions are straightforward. 


4. TIME DISCRETIZATION 

The resulting system of equations ( [T^ not only consist of ordinary differential equations (ODE), 
but also algebraic equations. The numerical time integration of a system of Differential Algebraic 
Equations (DAE) is more complex than the solution of an ODE. Therefore, we utilize the Newmark 
constant average time integration scheme with a special attention to select a consistent initial 
condition (Section|4T]). Based on the Newmark method, governing equations of ( fT^ is considered 
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at t. 


n+1 


Miln+i 

'^n +1 


Q '^n+l 


+ Mfgn +1 +-^'*^ 11+1 

ng 


Q Pn+\ — 1 


n+1 


■^^n+1 “t“ -^Qn+l ^ Pn+l —lF'n+1 

-^h 


(14) 


B^qa+i 


= 0 


in which the following assumptions are made to approximate the values at each time step 


Af . 

■Itn+l — U-a -\ -^Mn+1 

where At = tn+i — tn and 


++1 — Mn + ^++1 


At . ~ 2 . , 

Mn = Mn + Mn = “(Mn + 


gn +1 = + + ^9n+l 


gn = -(+ + ^gn) 


(15) 


using approximations of ( fTS] ! in the discretized equations of ([T^ results in the final time discretized 
version 

M Mn+l +Mf - Q = Pn+l 

++1 +^g„+l --BPn+1 ='Fn+i (16) 

Un+i +B^ g„+i = 0 

where the following definitions are made 


At^ 

M = M + ■ 


j ng At. . 


Q = Q 
At 


At 

T ’ 


B = B 


At 


(17) 


Pn+l = -Y (Pn+l - ^^itn - Mf Qn “ KUn) 
Pn+l = ^ (Pn+l - M+iin “ 


Rearranging and reverting Equation ( fT^ to matrix form yields 

■ M Mf -Q 
MfJ A -B 
-Q^ -B^ 0 



++1 


Pn+r 


gn +1 

= 

Pn+l 


Pn+l. 


0 


(18) 


0 represent a zero matrix of appropriate dimension. Above system of equations is symmetric with 
diagonal blocks M and A both symmetric positive definite. 


4.1. Consistent initial condition 


Following the three-field spatial approximation of the problem, the incompressibility condition, 
which is embedded in the mass conservation equation ([ 1 ^ 3 , forms an algebraic constraint. The 
presence of this constraint, makes ( fT^ represent a DAE system in time rather than an ODE in 
time. This fact is also emphasized by Diebels and Ehlers p7| . We note the absence of the algebraic 
variable (pressure) from the constraint ©a. We also note that not only the same matrix block 


-Q 

-B 


appears in (di. 2 and ([T 8 ]) 3 , but also the Schur complement of the upper left block matrix 


-[-Q^ 


■ M 

Mf 

-1 

-Q 


A 


-B 


is a nonsingular matrix with a bounded inverse. As a result, ( [T8| l yields to a Hessenberg index- 
2 system of DAEs as indicated by Ascher and Petzold |[72| and has the form of an augmented 
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Lagrangian equation. Consistency of initial conditions is one of the key aspects in DAE^ the initial 
values for variables in a DAE must satisfy not only the original equations in the system but also 
their differentials with respect to time | [74][75) . Accordingly, to select consistent initial conditions 
for acceleration terms and pressure, we must satisfy the original equations in the system as well as 
the differential of the mass conservation equation with respect to time, at time zero which reads 

Milo +Mfqo -Qpo=Po 

Mf^ilo + Aqo -BpQ=¥o (19) 

Q^ilo +B^qo =0 

Solving this equation, uo,qo,po is obtained. Note that in deriving [T9] fluid and solid phases are 
assumed to be initially at rest (-ito = 0,qo = 0) without any deformation (mq = 0 ), which is an 
intuitive set of initial conditions that also satisfy the mass balance equation at time zero. We observe 
that an inconsistent set of initial conditions results in instability of pressure and acceleration time 
histories illustrated by nonphysical oscillations as numerical artifacts. 


4.2. Equation of energy 

In this section, we derive the work done by various forces acting on the system during a time 
step. We show that the work done by pressure (Eagrange multiplier) vanishes, as expected for 
incompressible media, so the total mechanical energy of the system is preserved. 

The total incremental energy of the system is obtained by multiplying the discretized equations 
of motion (fT^i .2 at the midpoint n + ^ (i.e., average of the equations at time n and n + 1) by their 
respective displacement increments, and then adding them together. The resulting increment of the 
total energy is 

+ 5E^ + 5E^ + 5E^ = 6E^^ ( 20 ) 

where SE^'^ = E^^^ — E^^. Also dE^\ SE^\ SE^, SE^, 6E^, and SE^" denote the increments of 
kinetic energy of skeleton, kinetic energy of fluid, damping energy, strain energy of skeleton, work 
done by constraint force, and input energy respectively. Each of these terms are defined as 


SE'^’‘ := - iijMun) 


SE^‘ := Qn+i - ujMiqa) + -(q^^,Aqn+i - qjAqn) 


■= +9n)^^(9n+l +9n) := {uJ^yKUn+l -ujKUn) 

SE^'^ := — ((Mn+l + 'n„)^(P„+i + P„) + (Qn+l + 9n)^(®'n+l + 

SE'^ := + qJ^iB) + {ujQ + qjB)) p„^ i 

( 21 ) 

Where we employ the relationships of ( fT5] l for the displacement increments of solid and fluid. The 
statement of energy balance ( |20l i says that the increment of the input energy must be equal to the 
increment of the kinetic energy, strain energy and damping forces energy. We note that the increment 
of the work done by constraint force is zero, SE'^ = 0, since the terms inside the parentheses in 
(ID 4 are the discretized constraint equations that are set to be zero. Moreover, since we satisfy the 
constraint at initial time, we get E^ = 0 as well. 


§ An alternative approach to de al w ith incompressible poroelasticity system of DAEs is to utilize the energy preserving 
method proposed by Bauchau |73| . His approach satisfies the algebraic constraint at the end point (if written in terms 
of displacement), while imposes the algebraic variable at the midpoint; hence, it does not require an initial condition on 
the algebraic variable. This method reduces to the constant average Newmark method over unconstrained problems and 
it preserve the energyof the system by eliminating of the work done by algebraic variable. The result of our Numerical 
Examples in Sectionl^using the constant average Newmark method and Bauchau’s method perfectly match. 
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Using ( [^ , the energy components at each time step is recovered as 

'^n+l “ ^n+1 “ (^n+l-^f 9n+l) + ^(^n+l^Qn+l) 

+ SE^ Ell = <iKu„,i Eh = E^ + 6E^^ 

The damping and the input energy are time history dependent. 

5. IMPLEMENTATION 

The algorithm for implementation of the proposed numerical method for dynamic analysis of porous 
media is summarized in Procedure Procedure j^gives an outline of the construction of coefficients 


Procedure 1 Algorithm to solve dynamic poroelasticity_ 

1: Eorm the coefficient matrices, M, Mf, A, K , Q, B using Procedure © 

2: Select a time increment At (for an accurate result use Equation ( |27| l) 

3: Calculate M, A, Q, B using Equation ( [T7| ) 

4: Initialize ito, mq, Qq to satisfy @3 

5: Initialize a consistent pg to satisfy ([T^i ,2 and (03 [> an example is given in Section |4~T| 

6: for n 0, (numTimInc — 1) do 

7: Given (^^rsc)n 5 (2£rsc)n+i5 (Pprsc)n and (pprsc)n+i 5 Calculate IPn, IPn+i 5 and Fn-pi using 0. 

8: Update Pn +1 and Fn+i using Equation ( fT7| ) 

9: Solve 0 for Mn+I,< 7 n+ 1 > and 

10: Calculate rtn+i (if required evaluate accelerations as well) using ( [15] ) 

11: Eor the energy balance analysis: 

Calculate and using ( [2^ 

Calculate SEh and SEh using ( j^ 

Calculate, E^ = E? + SE^ and E^, = E^ + SE^ 

12 : end for 


defined in ( [T3] ) for a two dimensional problem. Section provides more details on the notations 
used here. JVfn and Afeg denote the number of free nodes and the number of free edges. In building 
these coefficient matrices, in addition to the connectivity (element to node) and nodal coordinates 
matrices used in the typical finite element problems, some additional data structures such as element 
to edge and node to edge are required. An example of the RT element implementation is elaborately 
discussed by Bahriawatiand and Carstensen ||76). To avoid clattering in Procedure]^ we introduce 


iVn,= 

VCm = 

Mm = 


Ni(x) 0 
0 JVi(x) 

1 


JVAx) 0 
0 hix) 


2A„ 


. i —3 or 6 for lineai* and quadratic triangle, respectively 
S^lll(x-Xi) Sm2l2{x- X2) - X3)] 


Nrn^N^dA , A^= WrhW^dA , iMtU= / N^^W^dA 


(23) 


These integrations over the element domains can be numerically evaluated by Gauss integration 
method.The following comments and Gauss quadrature rules over triangles |771 are considered for 
numerical examples of section]^ 


• The order of the Gauss integration must be selected to accurately calculate the integrations. 
Eor instance, at least a six-point formula (three-point formula), which gives accurate results 
for up to fourth (second) degree of precision, is required to approximate a mass matrix Mm 
with quadratic (linear) interpolation functions. 


Copyright © 2014 John Wiley & Sons, Ltd. 
Prepared using nmeauth.cis 


Int. J. Numer. Meth. Engng (2014) 
DOI: 10.1002/nme 










DYNAMIC POROELASTICITY 


13 


• A reduction in the order of Gauss numerical integration from the order that is required 
to evaluate these matrices exactly, leads to inaccurate results with long lasting and high 
amplitude nonphysical oscillations. 

• To decrease the computational effort in case of a quadratic interpolation function, we use 
a sub-parametric element rather than an iso-parametric one, since all the elements in the 
numerical examples are constructed as straight edges, and the geometry can be interpreted 
to a lower degree than the nodal variable’s interpolation. 

Moreover, in forming matrices of Procedure]^ we utilize the following properties of the RTq mixed 
element 

1 along edge j 
0 otherwise 

The right hand side equality is used in construction of matrix B. The left hand side equality indicates 
that the velocity interpolation functions at edge j, have the property of producing a unit flux through 
the j* edge and 0 flux through all other edges. Furthermore, according to our notation, in a typical 
triangle, the enumeration on the three vertex nodes is counterclockwise, while the edge enumeration 
is set such that the edge opposite to vertex node i is the edge of the triangle. 



v-iT, = 

0 otherwise 


(24) 


Procedure 2 Formation of coefficient matrices for two dimensional problems 
1 : Given definitions in ( fT3| ) 

2 : A •(— OATfeg ; M and K ^ Oatj,, ; Mf ^ Q ^ B ^ 

3: for m ^ 1, iVei do 

4: Assemble A^, and from ([2^ into A, M and Mf respectively 

5: Assemble Qm = [l 1 oj dA into Q 

6: Assemble ^ [s^i If Sm 2 k Sm 3 k] ^ dA into B 

7 : Assemble Km = X 4 Bm^CBmdA into K 

8 : end for 


6 . NUMERICAL EXAMPLES 


To demonstrate the capabilities of our formulation to simulate the incompressible dynamical 
behavior of saturated porous media, we perform two benchmark numerical examples and verify 
them with analytical solutions (Example 1) and boundary element solutions (Example 2). As was 


noted in section 3.2 it has been found that when the hydraulic conductivity ATh is very small, locking 


(in terms of failing the EBB condition) occurs To explore this effect for our formulation, 

we consider a third set of examples to evaluate the stability of the method when locking is often a 
concern. Overall, we found that the presented numerical scheme is stable and sufficiently accurate 
for a wide range of material parameters including large and small values of hydraulic conductivity. 


6.1. Example 1: half space — soil column 

The first benchmark example is the application of a uniform load on the surface of a half space 
|[^|5 ^57 ^78|. Due to symmetry, only a representative column of small width, 0.1 m, is 
studied and the problem is of one dimension. The model is taken as an incompressible saturated 
soil with plane strain behavior. The material data are shown in Table |I] and are taken from the 
same benchmark example as in |47 The geometry and finite elements are illustrated in Eigure 
1 ^ where the side walls and the bottom are impermeable and the displacements normal to their 
surfaces are constrained. The upper boundary is perfectly drained and subjected to a step load of 
3 ^ amplirnde, which results in a total load of 0.3kN on top boundary. 400 finite elements with a 
crisscross mesh pattern are used to model the soil column. We also use a time step of At = 1 x 10“"^ 


(see section 6.1.2i that is adequately small for the purpose of wave propagation analysis. We have 
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E 

V 

Ps 1^3 

kg 

Pf 5# 

n 


Example 1 

14.516e-i-3 

0.3 

2000 

1000 

0.33 

le-2 

Example 2 

14.516e-i-3 

0.3 

2700 

1000 

0.42 

le-1 and le-4 

Example 3 

10 

0.4 

2667 

1000 

0.4 

le-7 


Table I. Material properties for numercial examples 



Figure 4. Example 1 — a column of soil under step load: geometry and finite element mesh. We have 
isolated a 10 m length by 0.1 m width for analysis. The figure is not in scale. Three dots show the 

continuation. 


provided the results utilizing the Pi — RTq element and a lumped mass matrix unless otherwise 
mentioned. 

With this example, we pursue many goals as listed below: 


1. Verification with analytical solutions 

• Short-term behavior: we compare the numerical results with the analytical solution of 


de Boer et al. 1791 for wave propagation in an infinitely long column of incompressible 


saturated soil. The influence of the total depth of the column is negligible | [47| ; therefore, 
for our numerical analysis, we are able to isolate a 10-m length and still compare the 
short term behavior of our system with de Boer’s results for an infinite long column. 

• Long-term behavior: we compare long term behavior with Terzaghi’s equation for 
consolidation. 

2. Wave propagation analysis: according to the well known derivation of Biot l[TJ, there are three 
types of waves that pass through a porous medium. Two dilatational waves (fast and slow) 
and one shear wave[j It is proved analytically that as a consequence of the incompressibility 
condition, the fast dilatational wave propagates with an infinite speed and practically, only the 


slow dilatational wave exists 179 801. Also, note that for this one-dimensional problem, there 
is no shear waves. 

3. Pressure build-up and diffusion: we are interested in studying the initial and long term 
behavior of the pressure field. 

4. Effect of instantaneous loading: the impact and shock nature of the step load creates a 
challenging condition for the incompressible porous medium. The effect of this condition 
on the response history of the medium is studied, and some remedies are suggested to avoid 
the associated numerical artifacts. 


^Accounting for apparent mass, there are two shear waves, one in fluid and one in solid skeleton 


0 


Copyright © 2014 John Wiley & Sons, Ltd. 
Prepared using nmeauth.cis 


Int. J. Numer. Meth. Engng (2014) 
DOI: 10.1002/nme 







































DYNAMIC POROELASTICITY 


15 



Figure 5. Example 1 — displacement time history. Good agreement between numerical and the analytical 
computations is observed for both the transient behavior given by de Boer et al. ( 22 ) and the long term 


behavior given by Terzaghi 25 


5. Effect of mass representation: it is suggested in the literature, that such numerical artifacts 
can be alleviated by appropriate representation of mass matrix. We explore different mass 
lumping methods and compare their results with those of the consistent mass matrix method. 

6. Energy balance analysis: in order to highlight the working of the presented numerical 
algorithm, we monitor the balance of the energy in the porous medium. In particular, by 
means of the energy balance equation, we can evaluate our time integration scheme and size 
of the time step. 

6.1.1. Verification with analytical solutions 

Shown in Eigurej^is the displacement response at the top boundary. Our numerical results match 
well with de Boer’s analytical solution in the transient region | [79| . We point out that de Boer’s 
solution is for an infinitely long column where there is no reflection of wave from the bottom 
boundary; however, we isolate a finite length for computations, and therefore our numerical plots 
separates from de Boer’s solution after a while. Moreover, to verify the long term behavior, we use 
Terzaghi’s equation for steady state consolidation. According to that equation, the total settlement 
at the top boundary |[8T| equation 2-3] can be simplified for such one dimensional problems as: 


6s = rriy f L (25) 

where / is the total load, H is the length of the soil column, and TOv is the the slope from a one 
dimensional test, 

ds 

TOv = — = A + 2G (26) 

da 

where a plain strain assumption is made. 


6.1.2. Wave propagation analysis 
wave speed and choice of time step 

The significance of the choice of an adequate time step is clear in a wave propagation problem. In 
order to represent the wave travel accurately, it is recommended that the time step is determined by 


the approximate Courant-Eriedrichs-Eewy (CEE) condition (equation 9.102, |231) 


At < 


Ax„ 


Co 


(27) 


This is the smallest time that is required for a dilatational wave to pass across any of the elements. 
Given the smallest element size, Aatmin, we can compute At knowing the wave speed Cq. An 
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Figure 6. Example 1— snapshots of wave propagation through the medium. The propagation speed reduces 

over the time as roughly calculate in \29). 


estimate of the slow dilatational wave speed is given by de Boer et al. | [79| as, 

Co = 


rriy 




(28) 


Using the parameters of Example 1 in the above relationship, we get Co = 85.1™. In the next 
section, we compute wave velocity using the plots of Figure and demonstrate correlation with 
this estimate. 

Setting the minimum mesh size. Ax = ^ m, to be the short-edge length for triangles, equation 
( |Z7l i results in a time step of At < 8.3 x 10“'*, so we apply a time step of Af = 1 x 10“^. 

Wave velocity analysis 

We are interested in observing the passage of the slow dilatational wave which is clear in Figure]^ 
recorded at some snapshots of time. The time span was selected such that no reflections of the wave 
affect the results. In this plot, we see the propagation of the wave induced by the load at z = 10 m 
and its propagation through the medium more and more as time progresses. Using this figure, the 
approximate wave velocity can be calculated as: Cq = ^. For instance, taking the first three plots 
we get: 


Co 

Co 

Co 


10-7.8 „„m 

0.0025 - 0 “ s 

10 — 4.6 „„m 

-= 72 — 

0.075 - 0 s 

10-1.5 m 

- = 56— 

0.150-0 s 


(29) 


Recall that an approximate value of 85 f was calculated earlier using de Boer’s equation p8| ); that 
is close to the early time value of our numerical speed in Equation ( |29| . However, our numerical 
results show a decreasing trend of wave speed as the wave travels through the media which is 
predicted by Schanz and Pryl l ISO) as well. One possible reason for this difference is the finite length 
of the column in our example compared to infinite length of de Boer analytical formulation. 

Time histories of skeleton velocity at different locations (z), shown in Figure 8b also reveals 
the passage of the wave through the medium. We observe that moment when the wave reaches to 
a certain depth, the velocity profile of that depth separates from the rest of the profiles with zero 
velocity. 


6.J.3. Pressure build-up and dijfusion 

Pore fluid pressure time histories computed at different depths is displayed in Figure]^ for short and 
long term. Upon the instantaneous application of the load, a pressure build-up occurs at all points of 
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(a) (b) 

Figure 7. Example 1— (a) Pressure time histories. Diffusion of pressure is apparent over the long time, (b) 
Zoom in of first moments. Pressure jumps upon loading at time zero. Also note the diffusion begins right 

after the passage of the stress wave. 



Figure 8. Example 1— (a) Solid skeleton velocity time histories for some representative depths. The trend, 
however, is the same at all locations with their velocities ultimately reach to zero, (b) Zoom in of first 
moments at some critical depths. The high frequency oscillations, byproducts of step loading, are more 
pronounced close to the top boundary, but they subside as time marches forward. 


the medium, indicating the propagation of the pressure at infinite speed ||69| (Figure [7b|. After the 
stress wave hits a depth, it induces a change in pressure which results in diffusion of the pressure 
back toward its initial state. 

6.1.4. Effect of instantaneous loading 

The time history of the skeleton velocity is shown in Figure Over time, as energy dissipates 
and diffusion occurs, the velocity tends to zero (Figure]^. The trend in this plot is the same at all 
depths, so we only display the response at a few representative locations. 

The step loading due to its abruptness creates a jump in pressure and inertia forces. This in turn 
excites some higher modes of vibration that survive only for a short time interval due to presence 
of physical damping in the form of fluid diffusion. These early high frequency oscillations are 
pronounced in the vicinity of the top boundary with rapid pressure change as shown in Figure 
We recognize that these high frequency oscillations are a numerical artifact as their frequencies 
become unbounded upon mesh refinement. We note that decreasing the time step does not remove 
these spurious oscillations either. We also found that these nonphysical oscillations are not affected 
significantly (or removed) by finite element type, as presented in Figure In these plots, we 
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Figure 9. Example 1 — a comparison of the solid skeleton velocity time histories. In the proximity of the 
loading boundary, the P 2 — RTq element exhibits higher amplitudes and frequency oscillations compared to 
the linear Pi — RTq element. After a while, both elements converge to identical results. 


compare velocity response of the medium using Pi — RTg and P 2 — RTq elements with the same 
number of elements. It is observed that the response using P 2 — RTg element has higher amplitude 
and frequency oscillations in comparison with Pi — RTq element solution. However, as we move 
away form the loading boundary, the responses of both finite elements are indistinguishable. 

To avoid these numerical oscillations, we explore different remedies - (1) Replacing the step load 
with a ramp load (even with a relatively short rise time) smooths out the initial oscillations (the 
results of this study are not shown here); (2) An alternative approach is to apply time integration 
methods containing algorithmic dissipation such as the Generalized a method | (8^ . We note that the 
results shown here are computed by the constant average Newmark scheme where no algorithmic 
dissipation is used and consequently, the total energy of the system ( [20| ) is preserved. This fact is 
illustrated in Figure |llb| as well; (3) While such modifications of time integration algorithm can 
be easily implemented, we instead search for remedies at the element formulation level. In | [8^ , it 
is suggested that such oscillations may be associated with representation of mass matrix. This is 
discussed in the next section 6.1.5 We also point out that these early numerical oscillations do not 
appear in displacement time histories since the integration process automatically filters out the high 
frequency components. 


6.1.5. Ejfects of mass representation 

In constructing the matrix M defined in ( [T3] ), we use the same shape functions for velocities and 
displacements, resulting in a consistent variational formulation as presented in ( [TT] l. Accordingly, M 
represents a consistent mass matrix, which is generally non-diagonal. An alternative to the consistent 
mass matrix is to employ a diagonal mass matrix based on direct lumping. In this method, the total 
mass of element is directly apportioned to nodal DOFs, ignoring any cross coupling. There are 
several ways of constructing a diagonally lumped mass matrix. Here, we test two of them, namely, 
mass lumping by nodal Lobatto quadrature |84| and the lumping method proposed by Hinton | |85] . 
We note that in this example the lumping methods are applied on the mass matrices associated with 
the block related to the skeleton displacement, i.e, M and M{ [i] while we keep using the same 
mass matrix for the fluid velocity block, i.e., A; since there is no notion of fluid displacement in our 
formulation and instead the primary variable filed is fluid velocity. The resulting weighting factor of 
each DOF in the element mass matrix is illustrated in Figure p^a| Note that over the linear triangle 


II The lumping process can be viewed as substituting a different interpolation function for velocity (as opposed 
to using the same interpolation function as used for displacement) to approximate the mass matrix as, = 
Jy NiAN]^ dU. Accordingly, to build the coupled mass matrix Mf we must use the same interpolation function Nj^ 
in order to tune it with the lumped mass ■ 
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J_ _ 2 _ 

3 5 7 



linear triangular element quadratic triangulai' 
element: Hinton method 


(a) 


0 



quadratic triangular 
element: Lobatto method 



Time (s) 



Time {$) 


(b) 


(c) 


Figure 10. Example 1 — (a) Weighting factor of each nodal DOF in the element lumped mass matrix. Note 
that over the linear triangle element. Pi — RTq, both lumping procedures result in identical nodal weights, 
(b) Comparison of lumped and consistent mass matrix responses for the Pi — RTg element. The consistent 
mass matrix invokes a greater amount of high frequency oscillations, especially upon the arrival of the wave 
front, (c) Unstable response of the quadratic triangle element, P 2 — RTq, when a Lobatto quadrature lumping 

method is used. 


which is used in the Pi — RTq element, both procedures give rise to the identical nodal weights. 
Shown in Figure [TOb] we provide the numerical results employing both the lumped and consistent 
mass matrices for the linear triangle element. Pi — RTq. The lumped mass matrix results are 
displayed in solid line while the consistent mass matrix results are in dashed line. We see that the 
consistent mass matrix creates oscillations with higher frequencies compared to the lumped mass 
matrix response. In other words, it is more noisy compared to lumped mass response especially 
upon the arrival of the wave front. Nevertheless, in both cases, noise is damped out after a short 
while. 

With the quadratic triangle element, P 2 — RTq, we see the same trend as for linear triangle element 
when the Hinton et al. method of lumping is applied. However, the use of the Lobatto quadrature 
lumping method gives rise to an unbounded and unstable solution as plotted in Figure 10c This 
numerical instability occurs regardless of the mesh size, the time step size, and the hydraulic 
conductivity. In fact, we observe that the Lobatto quadrature lumping method over the quadratic 
triangular element produces unstable numerical solution even in a typical dynamic elasticity 
problem. We also emphasize that when using a quadratic triangle element, a consistent load vector as 
defined in ( fT3l l should be exerted regardless of mass matrix begin lumped or consistent. Otherwise, 
very inaccurate results are obtained as reported in as well. 


6.1.6. Energy plots 

We have developed the energy balance equation for the entire porous medium in section [T2| As a 
means of highlighting the performance and accuracy of the presented time integration approach, the 
relative error in the energy balance of the solution is computed as shown in Figure pTa| Left hand 
side (LHS) energy is defined as the summation of all the energy terms except the input energy as 
expressed in ( |2^ . The relative error is defined as the absolute value of the difference of the LHS 
energy and input energy divided by the input energy. The error is very small and the time step 
scheme is accurate. Also all the components of energy are illustrated in Figure [TTb| 
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(a) 


(b) 


Figure 11. Example 1 — (a) Energy balance error, (b) Components of energy as defined in l|22|> 


6.1.7. Summary of findings 

• As shown in Figure a good agreement between the numerical and the analytical results is 
observed for both the transient behavior given by de Boer et al. | |79) for wave propagation 
analysis, and the long term behavior given by Terzaghi’s consolidation equation. 

• We observe a decreasing trend of wave speed as it travels through the medium. 

• Pressure diffusion begins right after the passage of the stress wave as shown in Figure [7] 

• Due to the instantaneous nature of the step loading, high frequency numerical modes are 
excited early in time. These are shown in the velocity time history of Figure We suggest 
using a small ramp or a time integration method with algorithmic dissipation to alleviate this 
problem. Also a mass lumping method is shown to be effective in reduction of these types of 
numerical artifacts. We note that these high frequencies are not apparent in displacement time 
history. 

• The use of the Lobatto quadrature mass lumping method gives rise to an unbounded and 
unstable solution as shown in Figure [T0c| 

• The small relative error in the energy balance of numerical solution expresses the satisfactory 
performance and accuracy of the presented time integration method. 


6.2. Example 2: soil block with a partially loaded surface 

A key benchmark for verification of numerical scheme is the analysis of a saturated soil block with 


a partially loaded surface as shown in Figure 12 Due to symmetry only half of the problem is 
modeled. All the boundaries except the unloaded area of the upper boundary are impermeable. The 
displacements normal to the surfaces of the side walls and the bottom are constrained. We adopt 
the material properties that are used by Li et al. Q and shown in Table |T] which are also similar to 
those used by Diebels and Ehlers m Two different values of hydraulic conductivity is considered 
to cover the typical and the small range of this parameter. The model deforms in plane strain under 
a vertical step load of 15 ^ amplitude. The time increment is fixed at the value of 5 x 10“^. Other 
time steps including Af = 1 x 10“^, 1 x 10“'* are also tested but they appear to have no significant 
influence on the response. 

With this example as well, we pursue many goals as listed below: 


1. Verification with boundary element solutions 

• Significance of finite element parameters: We show that the dynamic behavior of porous 
media with small value of hydraulic conductivity is very sensitive to the choice of finite 
element, mesh pattern, and mesh size. 

2. Damping and frequency content: the dynamic response of the medium with different material 
parameters is studied in terms of damping and frequency content. 
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ISkPa 


node 4 
0 



side and bottom 
boundaries are 
impermeable 
(in blue) 


Figure 12. Example 2 — geometry and finite element mesh 


3. Fluid and pressure diffusion: to check the smoothness of the variable field results, the passage 
of flow through the porous medium and the pressure contours are illustrated at different 
moments of time. 


6.2.1. Verification with boundary element solutions 

Shown in Figure 13 are the vertical displacement time histories at node 1 and 2 (the left and 
right upper comer nodes shown in Figure \V2\ . Our numerical results with a crisscross pattern of 
meshing (shown in black color) match exactly with solutions obtained using an implementation 
of the boundary element formulation of Dargush and coworkers |57 The boundary element 
solutions are not shown in the figure because they overlap our solutions. We remark that several 
finite element parameters such as mesh pattern, mesh size, and type of element can significantly 
affect numerical results in terms of dissipation and frequency content. This topic is discussed in the 
next section l6!2.2l 


6.2.2. Significance of finite element parameters 

In order to demonstrate the effect of different parameters on the dynamic response of porous media, 
we analyze several types of mesh patterns (see Figure and sizes over our two finite elements, 
i.e.. Pi — RTq and P 2 — RTq. The results are summarized in Table|I^and are shown in Figure [T3lfor 
the case of Pi — RTq element. An important observation is made in Figure [T3a| for porous media 
with small hydraulic conductivity value - when using a criss pattern for griding, a coarse mesh 
would cause a dramatic numerical error in terms of dissipation and frequency content (shown in red 
color). However, after sufficiently refining the mesh (shown in blue color), the numerical solution 
converges to the correct one (boundary element solution). We remark that most of the studies, such 
as ||5 47 ^ present results similar to those shown in red color possibly due to insufficient mesh 
refinement and using low accuracy methods. 

On the other hand, in case of a typical value of hydraulic conductivity K^ = 1 x 10“^™, 
we observe that the dynamic behavior of the porous media is not significantly affected by 
aforementioned finite element parameters - only combination of a coarse mesh with criss pattern 
could slightly contaminate the frequency content and the amount of settlement (see red color results 
in Figure [T3b| which is expected due to an asymmetric pattern of meshing. In contrast to numerical 
results of the Pi — RTq element, the displacement time history results using the quadratic element 
P 2 — RTq, are satisfactory for all types of mesh patterns and sizes regardless of the hydraulic 
conductivity value. 

; the black color results using the crisscross pattern perfectly match the boundary element solution 
that is not shown here 
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element type mesh pattern mesh size performance Remarks 


criss 

Pi - RTo 

coarse 

X 

ATh = 1 X exhibit nonphysical 

dissipation and inaccurate frequency see 
Figure 13a iCh = lxlO“^“: slightly 
contaminate the frequency content and the 
amount of settlement 


fine 

/ 


crisscross & union jack 

coarse 

fine 

/ 

/ 



criss 

P 2 - RTo 

coarse 

/ 


fine 

/ 

In all four analysis, time history of veloc¬ 
ity develop high frequency oscillations 

criss ross & union jack 

coarse 

/ 

while displacement time history is smooth 

fine 

/ 


Table II. Example 2 — performance of the method using different element types, mesh patterns and mesh 

sizes 


- node 1 criss/coarse mesh 

- node 2 criss/coarse mesh 

- node 1 criss/fine mesh 

- node 2 criss/fine mesh 

-node 1 criss cross/coarse mesh 

-node 2 criss cross/coarse mesh 



(a) 



(b) 


Figure 13. Example 2 — displacement time histories for two different hydraulic conductivities, (a) Low 
hydraulic conductivity = 10“"^ numerical results are sensitive to finite element parameters such as 
mesh pattern and size. For example criss pattern with coarse mesh leads to erroneous results in terms of 
dissipation and frequency content (in red). Most of literature studies report such results. Mesh refinement 
fixes the problem (in blue), (b) Normal hydraulic conductivity Jfh = 10“^ numerical results are not very 
sensitive to finite element parameters; criss pattern with a coarse mesh did slightly skews the settlement and 
changes the frequency, whereas crisscross pattern works properly under a coarse mesh. 


6.2.3. Damping and frequency content — effects of hydraulic conductivity 

Comparing Figures 13a and |13b it is seen that the frequency of the dynamic response of the 
porous medium remains almost the same for models with different values of hydraulic conductivity, 
suggesting that the wave speed is not affected by this parameter. In terms of damping, when 
hydraulic conductivity is larger, the response decays faster in time and the steady state consolidation 
solution is achieved over a short period of time (Figure |13b| l due to the greater drainage and fluid 
diffusion. In contrast, in case of small hydraulic conductivity, oscillations at the corner nodes persist 
over a long period of time and the steady state consolidation solution is not achieved in that period. 
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(a) (b) 

Figure 14. Example 2 — pressure time histories, (a) — 1 x (b) K'a — 1 x 10“^ ”. 


In this condition, a nearly undrained response (similar to a nonporous incompressible elastic body) 
is observed where both the corner nodes move in skew symmetric manner. Also such low amount of 


drainage does not allow for any damping as shown in Figure 13a This behavior can be explained by 
considering the damping energy relationship of Equation (|21|i2. This energy is proportional not 
only to the damping coefficient but also to the fluid Darcy velocity. In case of small hydraulic 
conductivity iCh, although the damping coefficient increases, the fluid Darcy velocity drops to 
such an extent that results in the overall reduction of damping energy. Comparing the slope of 
the damping energy in Figures [TSb] and [T8b| reveals the same fact. 

6.2.4. Fluid and pressure diffusion 

Time history of the pressure is plotted in Figure [T4| for all the comer points. For a typical value of 


conductivity iTh = 1 x 10 ^ Figure 14b shows a general trend of pore pressure decay. The initial 


pressure build-up is apparent in this plot. On the other hand, the pressure time history for the small 
value of conductivity. Figure |14a[ shows a small amount of decay while oscillating with the same 
frequency as that of skeleton and fluid velocity. Shown in Figures and are several snapshots at 
different moments of time of the pore pressure contours and qualitative fluid velocity vector fields. 
The velocities in these figures are computed at the centroid of each triangular element using the 
interpolation function in ( [T0| ). 

Figure [^illustrates a comparison of fluid Darcy velocity and skeleton velocity. We observe that 
in an incompressible porous medium only the slow dilatational wave exists, which is characterized 
by motion of the fluid and skeleton with opposite phase angles | |80) . 

Furthermore, similar to Example 1, the satisfactory performance and accuracy of the time 
discretization technique is evaluated through the error in satisfying the energy balance of the system 
as shown in Figure [T^ 

6.2.5. Summary of findings 

• Numerical results are verified with boundary element solutions for a wide range of values of 
hydraulic conductivity (Figure [T3|). 

• The numerical simulation of dynamics of porous media with small value of hydraulic 
conductivity can be significantly biased with with an unsuitable choice of finite element 
method, size, or mesh pattern (Figure [T3a| . Most studies in the literature, such as ||^47 78) , 
present erroneous results in terms of dissipation and frequency content possibly due to 
insufficient mesh refinement or low accuracy methods. These results in the literature suggest 
an increase in damping with decreasing hydraulic conductivity, while a proper finite element 
discretization leads to the opposite result. 

• The frequency of the response remains almost the same for different hydraulic conductivity 
values (Figure [T^, suggesting that the wave speed is not affected by the value of the hydraulic 
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(a) 


(b) 


Figure 15. Example 2 — pressure contours and qualitative velocity vector field at A'h = 1 x 10“"^ The 
results are smooth. Velocity field is larger close to the common loading and free boundary on top. (a) t=0.1 

s. (b) t=0.9 s 




Figure 16. Example 2 — pressure contours and qualitative velocity vector field at ATh = 1 x 10 ^ (a) 

t=0.2 s. (b) t=l s. (c) t=2.2 s. (d) t= 3.8 s. 


conductivity. However, damping is shown to be dependent to the hydraulic conductivity value, 
i.e., the larger this value is, the response is damped out faster. In contrast, in case of small 
hydraulic conductivity (Figure [T3a| oscillations persist over a long period of time and a nearly 
undrained response (similar to a nonporous incompressible elastic body) is observed. 

• The slow dilatational wave in an incompressible porous medium is characterized by motion 
of the fluid and skeleton with opposite phase angles as shown in Figure [TT^ 
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Figure 17. Example 2 — skeleton and fluid Darcy velocity time histories, (a) = 1 x 10 (b) K\t = 

1 X 10“^ In an incompressible porous medium, the slow dilatational wave is characterized as moving of 

the fluid and solid in opposite phases. 


6.3. Example 3: challenge of locking - model with very low hydraulic conductivity 


As discussed in section 3.2 there are two limiting cases at which the finite element must satisfy 
the stability criterion. The rigid skeleton limit requires that flow variables (m,p) be compatible with 
each other in terms of LBB condition; Raviart-Thomas mixed finite element is theoretically proved 
to satisfy this condition. The other limiting problem is the very low hydraulic conductivity case 
(iCh —> 0) in which the system degenerates to an incompressible elastic system with u and p as the 
primary variables. For this limiting case, the finite element discretization of {u, p) must comply with 
LBB criterion. In this example we aim at studying such a limiting case under a special configuration 
that typically induces spurious pressure mode and locking. In this work, the phenomena of locking 
refers to the situation when a finite element discretization is not able to pass the inf-sup condition 
(see discussions in section |63T]|. 

This example seeks to fulfill the following goals 


1. Incompressible elasticity problem: we first examine the stability of finite element spaces of 
{u,p) in an incompressible elasticity problem {Kh = 0). We consider two elements, namely, 
the linear displacement-constant pressure Pi — Pq and the quadratic displacement-constant 
pressure P 2 — Po, and study each one on three mesh patterns as shown in Figure [T^ 

• Spurious modes and inf-sup test: the performances of various finite elements are 
explored under different scenarios: element level (local), assemblage level (with 
boundary conditions applied), and the uniform stability (assessed with increasing mesh 
refinement). 

2. Poroelasticity problem under very low hydraulic conductivity (iTh very small) 

• Checkerboard pattern in the pressure field: in a porous medium with very low hydraulic 
conductivity, numerical methods often produce checkerboard patterns and nonphysical 
oscillations in the pressure field. We examine the performance of our scheme over the 
time to address this issue from a dynamic point of view. 


6.3.1. Incompressible elasticity problem 

The full discretized matrix form of the incompressible elasticity problem can be extracted from 
by setting q = 0 as 


■ M -Q 


■*^11+1 


>■ 

_-QT 0 


Pn+1. 


0 


(HI' 

(30) 


Spurious modes and inf-sup test 

We investigate on the stability of our finite elements from three perspective: 
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Time(s} 

(a) 


- input energy 

-damping energy 

- strain energy 

-kinetic energy solid 

kinetic energy fluid 


0.25 




-0.05'- ^^^^-' 

0 2 4 6 8 10 

Time(s) 

(b) 


0.7 

0.6 

0.5 

0.4 

0.3 

0.2 

0.1 

0 


\J\fV^^ - 




|0.12 

- LHS energy 

. input energy 

- relative error 


-0.08 


-0.04i 

Q 

- 0.02 

-0 


2 4 6 8 10 

Time(s) 


-input energy 
-damping energy 

- strain energy 

- kinetic energy solid 
kinetic energy fluid 



(C) 


(d) 


Figure 18. Example 2 — energy plots. (a,b) Ky^ = lxlQ ™. (c,d) Ky^ = lxlQ ^ ™ 


A —1 m 


• Linear independence of constraint equation (with boundary conditions applied): one indicator 
of instability is the presence of spurious pressure modes, which occurs when the kernel of the 
Q is nontrivial. In this case, the spurious non-zero pressure mode satisfies 

QPs = 0 

This is the result of constraint equations being linearly dependent. Clearly, in the presence 
of these spurious modes, the inf-sup expression yields exactly zero and the problem is 
not solvable. The existence of such non-zero Ps can be crucially dependent on boundary 
conditions. 

• Linear independence of constraint equation (without boundary conditions applied): there is 
another type of spurious mode, which is the side effect of using certain mesh patterns. These 
are element-wise (local) spurious modes due to redundant constraints at the macroelement 
level regardless of any imposed boundary conditions. An extensive discussion on this issue is 
found in | |84) . 

We assess combination of several elements and mesh patterns in a two-dimensional plane 
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element type* 

mesh pattern^ 

A: element-wise 
linear independent^ 
(no BCs) 

B: linear 
independent^ 
(N=l) 

C: inf-sup 
condition^ 

(N —)■ 00 ) 


criss 

/ 

/ 

X 

1 

0 

crisscross 

X 

X 

X 


union jack 

/ 

X 

X 


criss 

/ 

/ 

/ 

to 

1 

0 

crisscross 

/ 

/ 

/ 


union jack 

/ 

/ 

/ 


Table III. Linear independency of the constraint in the element and global level (assemblage) as well as 
the numerical inf-sup test. The Pi — Pg crisscross macroelement contains element-wise and global spurious 
modes. Pi — Pg unio^ack macroelement has global spurious mode. 

* See Figure]^ ^ See Figure [l^-c, ^ See Figure [l^ 


strain problem and a cantilevered square block geometry (Figure 


1911. Table III shows the 


results of linear independence of the constraint block in equation i under various scenarios. 
Columns A and B report if the matrix of constraints, Q, has full column-rank. It is observed 
that the Pi — Pg crisscross macroelement has one redundant constraint in each macroelement, 
which causes an element-wise spurious pressure; this was first discovered by Mercier | |87) . 
We also note that the criss element does not contain any spurious pressure mode either in 
element level or in the assemblage level (when boundary conditions is applied). 

Uniform stability - inf-sup expression: in order to fully verify the stability condition, the inf- 
sup expression must be bounded from below by a constant value independent of mesh size. 
The inf-sup condition is not simply equivalent to linear independence of constraint equations 
which is a purely algebraic property, it is rather an analytic one that ensure the uniform linear 
independency as the mesh is increasingly refined f/O) . 

To perform the inf-sup test of Chapelle and Bathe | |8^ on our elements, we construct the 
matrix form of the inf-sup condition and normalize it. The corresponding eigenvalue problem 
is then solved and the smallest eigenvalue is considered as the inf-sup expression. We observe 
that the normalization by the displacement norm is inevitable otherwise the smallest singular 
values of the constraint matrix alone may never approach to any limit (as occurred for the 
P 2 — Pg element without normalization). It is noted that in the presence of spurious modes 
the inf-sup expression is exactly zero, however, we use the smallest non-zero eigenvalue to 
build the graph of Figure]^ We display in this figure the inf-sup expression evaluated for a 
series of five mesh refinements, N = 1,2,4,8,16 where N is the number of macroelements 
along each side (an example of iV = 4 is shown in Figure [T9|i). We observe that the Pi — Pg 
finite element violates the inf-sup condition no matter what mesh pattern is used as indicated 
by inf-sup values unbounded when N increases. In contrast, the P 2 — Pg element satisfies 
the inf-sup condition, indicated by asymptotically approaching of the inf-sup expression to a 
value greater than zero. Column C of Table [nl| summarizes the results of these tests. 


6.3.2. Poroelasticity under very low permeability 

By passing the inf-sup test, the P 2 — Pg element is shown to be capable of overcoming the 
problem of locking in the incompressible elasticity problem. We wish to establish that consequently, 
it performs well in the low-permeability poroelasticity problems as well. We note that, in 
poroelasticity the Pi — RTg (P 2 — RTg) element is the three field extension of Pi — Pg (P 2 — Pg) 
elements introduced earlier for incompressible elasticity. In this section we study the configuration 
of Figure 21 which is a 1 (m) x 1 (m) cantilever square bracket similar to the one given in the 
previous section for inf-sup analysis. This problem was first used by Liu pT) and later by Phillips 
and Wheeler | [5^ to demonstrate the problem of locking in poroelasticity. Plane strain behavior is 
assumed and the material properties, similar to those used by Liu p7|, are shown in Table To 
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Figure 19. Mesh patterns and the problem considered for the inf-sup test, (a) Criss pattern, (b) Crisscross 
pattern, (c) Union jack pattern, (d) Cantilever square bracket problem used for the inf-sup test, an instance 

of = 4 with a crisscross pattern is shown here. 



Figure 20. Numerical inf-sup test on finite elements. All P 2 — Pq elements satisfies the inf-sup condition, 

while the Pi — Pq elements do not. 



Figure 21. Example 3 — cantilever square bracket: Geometry and finite element mesh 


simulate the incompressible elasticity situation and induce locking, we study the case of very low 
permeability, iCh = 1 x 10“^. A step load of 1 ^ is applied uniformly on top and all the boundaries 
are impermeable. Also on the fixed left boundary, the normal and tangential displacements to the 
surface are constrained. 

Checkerboard pattern in the pressure field 

In Figures 22a and |22b| we present pressure distributions along cross section A-A (X = 0.25 m) 
for Pi — RTg and P 2 — RTq elements respectively. These are plotted at different instances of time. 
The nonphysical oscillations and checkerboard pattern in the pressure field in Figure [22a|reveals the 
problem of the Pi — RTq element in handling low hydraulic conductivity as also reported in ]37|56| . 
However, implementing a quadratic skeleton displacement field as in the P 2 — RTg element resolves 
the problem and results in a smooth pressure field (Figure [22b|. The pressure field at other cross 
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(a) Pressure along cross section A-A plotted at different 
instances of time using the Pi — RTq element 


(b) Pressure along cross section A-A plotted at different 
instances of time using the P 2 — RTq element 


X 10 ^ 




(c) Pressure along various cross sections 
at time=2.5 s using the Pi — RTq element 


(d) Pressure along various cross sections 
at time=2.5 s using the P 2 — RTq element 


Figure 22. Example 3 — nonphysical pressure oscillation is visible at different moments of time and cross 
sections when using the Pi — RTq element, whereas a smooth pressure field is obtained when using the 

P 2 — RTq element 


sections is plotted in Figure p2c] and 22d at time (t = 2.5), showing the same phenomena. To get 
further insight, the pressure over the whole domain is plotted in Figure]^ at different instances of 
time confirming similar observations. One indication of instability is the checkerboard patterns on 
the pressure field with the Pi — RTq element (left figures), whereas the P 2 — RTq element (right 
figures), which satisfies the inf-sup condition, effectively eliminates locking. 

Dynamic analysis affects oscillations: 

We observe in Figure [22a| the incapability of Pi — RTq element in computing a correct pressure field, 
with the checkerboard patter persisting over time. In contrast, in quasistatic poroelasticity, it has 
been observed that these instabilities are pronounced in the beginning, but tend to disappear as time 
progresses p8]|56||6^. Such behavior in dynamic analysis with low hydraulic conductivity could be 


explained by the periodic deformation and velocity time histories; in such periodic condition, one 
encounter the static problem over and over, leading to a frequent instability issue. One a different 
note, by comparing Figure [2^ with [2^ we see that in early moments of time (up to about f = 0.1 
s) both Pi — RTq and P 2 — RTq elements predict similar and smooth pressure fields. One possible 
reason could be the sudden application of step load, which initially creates a dominant skeleton 
kinetic energy as illustrated in Figure]^ 


6.3.3. Summary of findings 

• The stability of our finite elements in terms of local and global spurious modes as well as 


the inf-sup condition is studied (Table IIIi. It is found that the crisscross mesh pattern on 
the Pi — Pq element contains a redundant constraint, which causes an element-wise spurious 
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X(m) 


X 10"^ 




6 


5 


4 


I 


(a) time=0.1 s 




X(m) 

(c) time=0.3 s 


X(m) 




3 

2 



(e) time=4.1 s 




(g) time=6.3 s 



(d) time=0.3 s 




(h) time=6.3 s 


Figure 23. Example 3 — the left pressure contours illustrate the results of the Pi — RTq element and its 
associated locking response with the checkerboard pattern, while the right ones are for the P 2 — RTg stable 
and locking free elements. Note that the colorbar is not identical in all plots. 
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Time(s) 

Figure 24. Example 3 — components of energy. Initially, the kinetic energy of s keleto n is dominant that 
prevents the locking even in unstable elements shown in Figure [2^ 


pressure. Also, results indicate that the Pi — Pq element violates the inf-sup condition, 
whereas the P 2 — Pq element satisfies it. 

The checkerboard pattern and nonphysical oscillation in the pressure field, shown in Figures 
|22a| and left Figures of 23 reveals the problem of the Pi — RTq element in handling dynamic 
analysis of porous medium with low hydraulic conductivity. Therefore, we introduce a 
quadratic skeleton displacement field in the P 2 — RTq element. This resolves the problem 
and results in a smooth pressure field (Figure [22b| and right Figures in[23]l. 

In dynamic analysis, the checkerboard patterns produced by an unstable finite element is 
shown to persist over time and is not eliminated, in contrast to static analysis, where these 
instabilities appear early but disappear at later time. 


7. CONCLUDING REMARKS 

An accurate and stable numerical framework for dynamics of incompressible saturated porous 
media is presented and verified over a variety of benchmark numerical studies. Our topology- 
motivated finite element is based on coupling of the mixed RTq and the nodal Galerkin elements, 
implemented based on a three-field {u — w — p) formulation. As a result, we do not neglect 
fluid acceleration terms, which enables applying the method to problems with considerable fluid 
dynamics such as soil liquefaction 0 and biomechanics of porous tissues under rapid external 
loading pO][T^ . We also demonstrate stability of the method regarding the EBB condition in the 
limiting cases (rigid skeleton and very low hydraulic conductivity). In particular, the P 2 — RTq 
element, which uses a quadratic Galerkin discretization for the skeleton displacement, overcomes 
locking, which otherwise manifests in the form of checkerboard patterns in the pressure field. We 
also observe that with an unstable finite element, checkerboard patterns can persist over time. Given 
that we directly calculate the physical damping in the form of fluid diffusion, our numerical results 
suggest a decrease in damping with decreasing hydraulic conductivity, while frequency content is 
not affected by hydraulic conductivity. It is noteworthy that an appropriate prediction of dynamic 
behavior of porous with small hydraulic conductivity requires specific attention to the choice of 
finite element, mesh pattern, and mesh size. 
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NOTATION 


H(div; O) 

H^n) 

L^{n) 

N{x), W{x 

iVfeg 

iVfn 


Kh, rif 

Pi - RTo 
P2 - RTn 


W,ix) 


iT(div; O) = {w : w e (T^(0))2, V • w G ^^(O)} 

= {u : u G Du G L^(0)} 

square integrable functions in O 

I, H{x) skeleton displacement, fluid velocity, and pressure interpolation functions 
number of free edges 
number of free nodes 
vector of pressure DOFs 
vector of Darcy velocity DOFs 
vector of skeleton displacement DOFs 

unit normal vector to edge j chosen with a global fixed orientation 
outward unit normal vector to edge j of element m 
velocity vector defined on edge j 
hydraulic conductivity, porosity 

linear Galerkin finite element coupled with RTq finite element 
area of element m 

quadratic Galerkin finite element coupled with RTq finite element 
average density of the medium, fluid density 
Lowest order Raviart-Thomas element 
effective stress 

velocity interpolation function for edge j 
prescribed traction 


u, w, p 


length of edge j 

prescribed pressure 

sign convention equals to hmj ■ Sj 

skeleton displacement, Darcy velocity, pressure 
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